function udt = h2eom(t, u, params)
%   choice of Pi: 
%   Pi_xx =  2 P1 (1+UX^2) +  2 UX UY P2
%   Pi_yy = -2 P1 (1+UY^2) +  2 UX UY P2
%   Pi_xy =  P2 ( 2 + UX^2 + UY^2 )

    A = params.A;
    Nx = params.Nx;
    Ny = params.Ny;
    Nt = Nx.*Ny;


    % unpack variables 
    T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
    ux = reshape(u(1*Nt+1 : 2*Nt), Nx, Ny);
    uy = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
    P1 = reshape(u(3*Nt+1 : 4*Nt), Nx, Ny);
    P2 = reshape(u(4*Nt+1 : 5*Nt), Nx, Ny);

    Tdx  = dx(T);
    Tdy  = dy(T);
    uxdx = dx(ux);
    uxdy = dy(ux);
    uydx = dx(uy);
    uydy = dy(uy);
    P1dx = dx(P1);
    P1dy = dy(P1);
    P2dx = dx(P2);
    P2dy = dy(P2);

    m11 = (-3).*T.^2.*(1+ux.^2+uy.^2).^(1/2);
    m12 = (1/2).*(1+ux.^2+uy.^2).^(-1/2).*((-3).*T.^3.*ux+(-4).*P1.*ux.*(1+uy.^2)+(-2).*P2.*uy.*(2+(-1).*ux.^2+uy.^2));
    m13 = (1/2).*(1+ux.^2+uy.^2).^(-1/2).*(((-3).*T.^3+4.*P1.*(1+ux.^2)).*uy+(-2).*P2.*ux.*(2+ux.^2+(-1).*uy.^2));

    m22 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(2.*P2.*ux.*uy.*(2+3.*ux.^2+3.*ux.^4+uy.^2+(-3).*uy.^4)+3.*T.^3.*(ux.^4+3.*ux.^2.*(1+uy.^2)+2.*(1+uy.^2).^2)+4.*P1.*(ux.^4.*(1+(-3).*uy.^2)+2.*(1+uy.^2)+(-3).*ux.^2.*((-1)+uy.^4)));
    m23 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(ux.*uy.*((-3).*T.^3.*(1+ux.^2+uy.^2)+4.*P1.*(1+ux.^2).*(1+3.*ux.^2+3.*uy.^2))+P2.*(8+(-6).*ux.^4+(-6).*ux.^6+12.*uy.^2+8.*uy.^4+ux.^2.*(8+6.*uy.^2+6.*uy.^4)));
    m24 = 2.*(ux+ux.^3).*(1+ux.^2+uy.^2).^(-1/2);
    m25 = uy.*(1+ux.^2+uy.^2).^(-1/2).*(2+3.*ux.^2+uy.^2);

    m32 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(2.*P2.*(4+4.*uy.^2+(-3).*uy.^4+(-3).*uy.^6+3.*ux.^2.*(2+uy.^2)+ux.^4.*(4+3.*uy.^2))+(-1).*ux.*uy.*(3.*T.^3.*(1+ux.^2+uy.^2)+4.*P1.*(1+uy.^2).*(1+3.*ux.^2+3.*uy.^2)));
    m33 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(2.*P2.*ux.*uy.*(2+ux.^2+(-3).*ux.^4+3.*uy.^2+3.*uy.^4)+3.*T.^3.*(2+2.*ux.^4+3.*uy.^2+uy.^4+ux.^2.*(4+3.*uy.^2))+4.*P1.*((-2)+(-3).*uy.^2+3.*ux.^4.*uy.^2+(-1).*uy.^4+ux.^2.*((-2)+3.*uy.^4)));
    m34 = (-2).*(1+ux.^2+uy.^2).^(-1/2).*(uy+uy.^3);
    m35 = ux.*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+3.*uy.^2);

    m42 = (1/4).*T.^2.*ux.*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+3.*uy.^2);
    m43 = (-1/4).*T.^2.*uy.*(1+ux.^2+uy.^2).^(-1/2).*(2+3.*ux.^2+uy.^2);

    m52 = (1/2).*T.^2.*uy.*(1+uy.^2).*(1+ux.^2+uy.^2).^(-1/2);
    m53 = (1/2).*T.^2.*ux.*(1+ux.^2).*(1+ux.^2+uy.^2).^(-1/2);


    R1 = (-3).*T.^2.*Tdx.*ux+(-3).*T.^2.*Tdy.*uy+(1/2).*uxdx.*(1+ux.^2+uy.^2).^(-1).*(2.*P2.*ux.*(ux+(-1).*uy).*uy.*(ux+uy)+(-4).*P1.*(1+ux.^2).*(1+uy.^2)+(-3).*T.^3.*(1+ux.^2+uy.^2))+(-1).*uxdy.*(1+ux.^2+uy.^2).^(-1).*(2.*P1.*ux.*uy.*(1+uy.^2)+P2.*(2+ux.^2+(-1).*((-3)+ux.^2).*uy.^2+uy.^4))+(1+ux.^2+uy.^2).^(-1).*(2.*P1.*ux.*(1+ux.^2).*uy+(-1).*P2.*(2+ux.^4+uy.^2+(-1).*ux.^2.*((-3)+uy.^2))).*uydx+(1/2).*(1+ux.^2+uy.^2).^(-1).*(4.*P1.*(1+ux.^2).*(1+uy.^2)+2.*P2.*ux.*uy.*((-1).*ux.^2+uy.^2)+(-3).*T.^3.*(1+ux.^2+uy.^2)).*uydy;
    R2 = (3/2).*T.^2.*Tdx+2.*P1dx.*(1+ux.^2)+2.*P2dx.*ux.*uy+P2dy.*(2+ux.^2+uy.^2)+(1/4).*uxdx.*(1+ux.^2+uy.^2).^(-1).*(3.*T.^3.*ux.*(1+ux.^2+uy.^2)+4.*P1.*ux.*(1+uy.^2+ux.^2.*(1+(-3).*uy.^2))+2.*P2.*uy.*(ux.^2.*(4+3.*ux.^2+(-3).*uy.^2)+4.*(1+uy.^2)))+(1/2).*uxdy.*(1+ux.^2+uy.^2).^(-1).*(P2.*ux.*((-2)+ux.^2+((-5)+3.*ux.^2).*uy.^2+(-3).*uy.^4)+3.*uy.*((-2).*P1.*ux.^2.*(1+uy.^2)+T.^3.*(1+ux.^2+uy.^2)))+(1/2).*ux.*(1+ux.^2+uy.^2).^(-1).*(6.*P1.*ux.*(1+ux.^2).*uy+P2.*((-2)+(-3).*ux.^4+uy.^2+ux.^2.*((-5)+3.*uy.^2))).*uydx+(1/4).*(1+ux.^2+uy.^2).^(-1).*(12.*P1.*(ux+ux.^3).*(1+uy.^2)+(-3).*T.^3.*ux.*(1+ux.^2+uy.^2)+2.*P2.*uy.*(4+4.*ux.^2+(-3).*ux.^4+(4+3.*ux.^2).*uy.^2)).*uydy;
    R3 = (3/2).*T.^2.*Tdy+2.*P2dy.*ux.*uy+(-2).*P1dy.*(1+uy.^2)+P2dx.*(2+ux.^2+uy.^2)+(1/2).*uxdy.*uy.*(1+ux.^2+uy.^2).^(-1).*((-6).*P1.*ux.*uy.*(1+uy.^2)+P2.*((-2)+ux.^2+((-5)+3.*ux.^2).*uy.^2+(-3).*uy.^4))+(1/4).*uxdx.*(1+ux.^2+uy.^2).^(-1).*((-3).*uy.*(4.*P1.*(1+ux.^2).*(1+uy.^2)+T.^3.*(1+ux.^2+uy.^2))+2.*P2.*ux.*(4+4.*uy.^2+(-3).*uy.^4+ux.^2.*(4+3.*uy.^2)))+(1/2).*(1+ux.^2+uy.^2).^(-1).*(3.*T.^3.*ux.*(1+ux.^2+uy.^2)+uy.*(6.*P1.*ux.*(1+ux.^2).*uy+P2.*((-2)+uy.^2+ux.^2.*((-5)+(-3).*ux.^2+3.*uy.^2)))).*uydx+(1/4).*(1+ux.^2+uy.^2).^(-1).*((-4).*P1.*(1+ux.^2).*uy+4.*P1.*((-1)+3.*ux.^2).*uy.^3+3.*T.^3.*uy.*(1+ux.^2+uy.^2)+2.*P2.*ux.*(4+4.*uy.^2+3.*uy.^4+ux.^2.*(4+(-3).*uy.^2))).*uydy;
    R4 = (1/2).*T.^2.*ux.*uxdy.*uy+P1.*(2+ux.^2+uy.^2)+(1/4).*T.^2.*uxdx.*(2+ux.^2+uy.^2)+(-1/2).*T.^2.*ux.*uy.*uydx+(-1/4).*T.^2.*(2+ux.^2+uy.^2).*uydy;
    R5 = (1/2).*T.^2.*uxdy.*(1+uy.^2)+P2.*(2+ux.^2+uy.^2)+(1/2).*T.^2.*(1+ux.^2).*uydx;

    mdet = m11.*(m25.*m34+(-1).*m24.*m35).*(m43.*m52+(-1).*m42.*m53);



    i11 = (m25.*m34+(-1).*m24.*m35).*(m43.*m52+(-1).*m42.*m53);
    i14 = (-1).*(m25.*m34+(-1).*m24.*m35).*(m13.*m52+(-1).*m12.*m53);
    i15 = (m25.*m34+(-1).*m24.*m35).*(m13.*m42+(-1).*m12.*m43);

    i24 = m11.*((-1).*m25.*m34+m24.*m35).*m53;
    i25 = m11.*(m25.*m34+(-1).*m24.*m35).*m43;

    i34 = m11.*(m25.*m34+(-1).*m24.*m35).*m52;
    i35 = m11.*((-1).*m25.*m34+m24.*m35).*m42;

    i42 = m11.*m35.*((-1).*m43.*m52+m42.*m53);
    i43 = m11.*m25.*(m43.*m52+(-1).*m42.*m53);
    i44 = m11.*((-1).*m25.*m33.*m52+m23.*m35.*m52+m25.*m32.*m53+(-1).*m22.*m35.*m53);
    i45 = m11.*(m25.*m33.*m42+(-1).*m23.*m35.*m42+(-1).*m25.*m32.*m43+m22.*m35.*m43);

    i52 = m11.*m34.*(m43.*m52+(-1).*m42.*m53);
    i53 = m11.*m24.*((-1).*m43.*m52+m42.*m53);
    i54 = m11.*(m24.*m33.*m52+(-1).*m23.*m34.*m52+(-1).*m24.*m32.*m53+m22.*m34.*m53);
    i55 = m11.*((-1).*m24.*m33.*m42+m23.*m34.*m42+m24.*m32.*m43+(-1).*m22.*m34.*m43);


    i11 = i11./mdet;
    i14 = i14./mdet;
    i15 = i15./mdet;

    i24 = i24./mdet;
    i25 = i25./mdet;

    i34 = i34./mdet;
    i35 = i35./mdet;

    i42 = i42./mdet;
    i43 = i43./mdet;
    i44 = i44./mdet;
    i45 = i45./mdet;

    i52 = i52./mdet;
    i53 = i53./mdet;
    i54 = i54./mdet;
    i55 = i55./mdet;

    
    Tdt  = -(i11.*R1+i14.*R4+i15.*R5);
    uxdt = -(i24.*R4+i25.*R5);
    uydt = -(i34.*R4+i35.*R5);
    P1dt = -(i42.*R2+i43.*R3+i44.*R4+i45.*R5);
    P2dt = -(i52.*R2+i53.*R3+i54.*R4+i55.*R5);


    udt = zeros(Nx.*Ny.*5, 1);

    udt(0*Nt+1 : 1*Nt) = Tdt(:);
    udt(1*Nt+1 : 2*Nt) = uxdt(:);
    udt(2*Nt+1 : 3*Nt) = uydt(:);
    udt(3*Nt+1 : 4*Nt) = P1dt(:);
    udt(4*Nt+1 : 5*Nt) = P2dt(:);


end
